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Exact simulation of max-stable processes 

Clement Dombry* Sebastian Engelke^ and Marco Oesting* 


Abstract 


Max-stable processes play an important role as models for spatial extreme events. Their 
complex structure as the pointwise maximum over an infinite number of random functions 
makes simulation highly nontrivial. Algorithms based on finite approximations that are used 
in practice are often not exact and computationally inefficient. We will present two algo¬ 
rithms for exact simulation of a max-stable process at a finite number of locations. The 


first algorithm generalizes the approach by Dieker and Mikosch [2015] for Brown-Resnick 


processes and it is based on simulation from the spectral measure. The second algorithm 
relies on the idea to simulate only the extremal functions, that is, those functions in the con¬ 
struction of a max-stable process that effectively contribute to the pointwise maximum. We 
study the complexity of both algorithms and prove that the second procedure is always more 
efficient. Moreover, we provide closed expressions for their implementation that cover the 
most popular models for max-stable processes and extreme value copulas. For simulation 
on dense grids, an adaptive design of the second algorithm is proposed. 


Keywords: exact simulation; extremal function; extreme value distribution; max-stable process; 
spectral measure. 


1 Introduction 


Max-stable processes have become widely used tools to model spatial extreme events. Occurring 
naturally in the context of extremes as limits of maxima of independent copies of stochastic 


processes, they have found many applications in environmental sciences; see for instance Coles 
fll993] , |Buishand et al.|P008| , |Blanchet and Davison| [ |20TT| , |Davison et al.|Q20 1 2Q . 

Any sample continuous max-stable process Z with unit Frechet margins on some compact do¬ 
main X C is characterized by a point process representation [de Haan 1984] 


Z(x) = max(t^(i), x G X, (1) 

i> 1 
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Figure 1: The Poisson point process > 1} (grey)- Only finitely many (colored) 

(£i, ipi) contribute the maximum process Z (bordered in black). 


where {((i, V’i), * > 1} is a Poisson point process on (0, oo) x C with intensity measure C _2 dC x 
v(<hji) for some locally finite measure v on the space C = C(X, [0, oo)) of continuous non¬ 
negative functions on X such that 


Jc^i x M dVO = 1, 


x e x. 


( 2 ) 


Figure [I] shows a realization of Z as a mixture of different random functions of the above point 
process. Due to this complex structure of max-stable processes, in many cases, analytical ex¬ 
pressions are only available for lower-dimensional distributions and related characteristics need 
to be assessed by simulations. Moreover, non-conditional simulation is an important part of con¬ 
ditional simulation procedures that can be used to predict extreme events given some additional 
information [see Dombry et al.[ 2013[ Oesting and Schlather 2014 for example]. Thus, there is 
a need for fast and accurate simulation algorithms. 

As the spectral representation (jTj) involves an infinite number of functions, exact simulation 
of Z is in general not straightforward and finite approximations are used in practice. For the 


widely used Brown-Resnick processes [Kabluchko et al. 20091, for instance, Engelke et al. 


[20111 and Oesting et al. [20121 exploit the fact that the representation (jT|) is not unique in order 
to propose simulation procedures based on equivalent representations. However, often these 
approximations do not provide satisfactory results in terms of accuracy or computational effort. 
Exact simulation procedures can so far be implemented only in special cases. Schlather] [20021 
proposes an algorithm that simulates the points {Q,i > 1} in (JTj) subsequently in a descending 
order until some stopping rule takes effect. If v is the probability measure of a stochastic process 
whose supremum on X is almost surely bounded or if Z is mixed moving maxima process with 
uniformly bounded and compactly supported shape function, this procedure allows for exact 
simulation of Z. For extremal-t processes [OpitzJ 20131, the elliptical structure of Gaussian pro¬ 
cesses can be exploited to obtain exact samples [ Thibaud and Qpitz]|2014[ . |Oesting et al.| [ |20 1 3] | 
focus on a class of equivalent representations for general max-stable processes that, in principle, 
allow for exact simulation in an optimal way in terms of efficiency. They propose to simulate 
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max-stable processes via the normalized spectral representation with all the spectral functions 
sharing the same supremum. Being efficient with respect to the number of spectral functions, 
the simulation of a single normalized spectral function might be rather intricate in some cases 
including the class of Brown-Resnick processes. For the latter, |Dieker and Mikosdi| [2015] 
recently proposed a representation that enables exact simulation at finitely many locations. 
Besides, several articles focus on the simulation of finite dimensional max-stable distributions 


or, equivalently, of their associated extreme value copula. Ghoudi et al. [1998] and Caperaa 


et al. [2000J propose simulation procedures for certain bivariate extreme value distributions. 


Stephenson [ |2003| considers the case of extreme value distributions of logistic type. Boldi 


[2009] provides a method for exact simulation from the spectral measure of extremal Dirichlet 
and logistic distributions. 

In this paper, we propose two methods to simulate a general max-stable process Z exactly at 
a finite number of locations. At first, we propose a generalization of the algorithm by Dieker 


and Mikosch 12015] showing that their approach relies on sampling from the spectral measure 


on the T|-sphere of a multivariate extreme value distribution. The main idea of the second 
procedure is to simulate out of the i nfinite set > 1} only the extremal functions [cf. 


Dombry and Eyi-Minko 2012 |2013|, i.e. those functions that satisfy Qwfx) = Z(x ) for some 


x G X (the colored functions in Fig. Q]). In contrast to all existing simulation procedures, the 
process Z is not simulated simultaneously, but subsequently at different locations, rejecting all 
those functions that are not compatible with the process at the locations simulated so far. Both 
new procedures are based on random functions following the same type of distribution that can 
be easily simulated for most of the popular max-stable models. Our algorithms also apply very 
efficiently to exact simulation of finite-dimensional max-stable distributions or, equivalently, of 
the associated extreme value copulas. 


2 Extremal functions 


Without loss of generality, we may henceforth assume that Z is a sample-continuous pro¬ 
cess with unit Frechet margins given by the spectral representation ([!]). Indeed any sample- 
continuous max-stable process can be obtained from a process with unit Frechet margins via 
marginal transformations. We provide in this section some preliminaries on extremal functions 
and their distributions that will be essential in the new simulation methods. We use a point pro¬ 
cess approach and recall first that the C-valued point process $ = {0?}?:> i with o l = Qipi is a 
Poisson point process with intensity 


d( A ) = fc fo° 1 fCV’6A}C 2 dC J'(dV’), dec Borel. (3) 


Definition 1. Let K e .1' be ci nonempty compact subset. . \ function (!) G $ is called ft —extremal 
if there is some x G K such that oix) = Z(x), otherwise it is called K-subextremal. We denote 
by the set of K-extremal functions and by $ f the set of K-sub extremal functions. 


It can be shown that and ( 1> K are properly defined point process. When K = {.z' 0 }, x 0 e X , 
is reduced to a single point, it is easy to show that ^ is also almost surely reduced to a single 
point that we denote by d>f o , termed the extremal function at point xo. The distribution of 0+ is 
given in the next proposition [see Proposition 4.2 in Dombry and Eyi-Minko, 20131. 
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Proposition 1. The random variables Z{x o) and /Z[x o) are independent. Furthermore, 

Z(x o) has a unit Frechet distribution and the distribution of 0t (] /Z( xq) is 

Pxa(A) = pr(0+/Z(x o ) G A) = J l {f/ f(xo)£A}f(x 0 ) u(d /), A c C Borel. (4) 

By definition, ff o (x 0 ) = Z(x 0 ). This entails that the distribution P Xo is supported by the subset 
of functions {/ G C, f(x 0 ) = 1}. 

Proposition 2. The restricted point process 

$ n {/ e c, f(xo) > 0} 


is a Poisson point process with intensity 


/ l{/(x 0 )>o}M<l/) = / / !{C/6A}C dC P XQ (df), A C C Borel. 

IA JC Jo 


(5) 


Proof. The fact that the restricted point process $ D {/ G C, f(x 0 ) > 0} is a Poisson point 
process with intensity l{/(x 0 )>o}Md/) is standard. We prove Equation ([5]). For A C C Borel, 



c Jo 


l {c/eA} C dC^ 0 (d/) = 



C Jo 


l {U/f{x 0 )eA}C ' dC f(x 0 ) u(df) 



IjCVcAjC dC l{/(s 0 )>0} u {df) — / l{/eA}l{/(a:o)>0}Af(d/). 

C Jo Jc 

Here, we use successively Eq. ([4]), the change of variable C = C/ f( x o) with f(x 0 ) > 0 and Eq. 
([3]) for the last equality. □ 

REMARK 1 . As a consequence of ©, independent copies Y lt i > 1, of processes with distri¬ 
bution P Xo result in a point process { C*Tj:}?;>i which has the same distribution as the restricted 
point process D {/ G C, f(x 0 ) > 0}. If u({f G C, f(x 0 ) = 0}) = 0, then <f> consists only of 
functions with positive value at Xq and <I> has the same distribution as {CiYi}i>i- This provides 
an alternative point process representation of the max-stable process Z in terms of a random 


process Y such that Y (xo) = 1 almost surely. In Engelke et al. 12014] and Engelke et al. [2015], 
this representation is exploited for statistical inference of Z. 

It follows clearly from Definition[T]that we have the decomposition 4> = The following 

proposition will play a crucial role in our second approach based on extremal functions. If fi, / 2 
are two functions on X, the notation /i < K f 2 means f (x) < / 2 (x) for all x G K. 


Proposition 3 (Dombry and Eyi-Minko [2012], Lemma 3.2). The conditional distribution of 
with respect to is equal to the distribution of a Poisson point process on C with intensity 

1{ f< K z}P{df )• 
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3 Exact simulation procedures 


3.1 Introduction 


Recently, Dieker and Mikosch [120151 provided a new representation of stationary Brown- 
Resnick processes that allows for exact simulation of their finite-dimensional distributions. In 
this section we introduce two methods for exact simulation of arbitrary max-stable processes 
and distributions. More precisely, for a fixed number N G N of pairwise distinct locations 
x — (xi, ..., xn ) G X n , we aim at obtaining exact simulation of the max-stable random vector 


Z(x) = (Z(x i),..., Z(x N )). 


( 6 ) 


Both procedures are intimately connected with the distribution P x in ([4]). 

The first method extends the Dieker and Mikosch ]2015[ approach. In fact, we show that their 
representation is nothing else than the so-called spectral representation of the Brown-Resnick 
process, and their procedure actually enables exact simulation from the spectral measure on the 
L|-sphere. In Section [32] we derive a way of simulating from the spectral measure of a general 
max-stable distribution or a possibly non-stationary max-stable process. 

The second procedure presented in Section [373] relies on conditional distributions of the Poisson 
point process underlying any max-stable process. This approach also allows for exact simulation 
of ([6]) by simulating at each location only the unique function that actually attains the maximum; 
see Fig. jT] It is intuitive and turns out to be even more powerful than the spectral method. 


3.2 Simulation via the spectral measure 

Let us recall the spectral decomposition of the max-stable random vector Z(x); for details we 
refer to Resnick [2008, Chap. 5]. Here, we write f(x) = (f(x i),..., /(xjv)) for the restriction 
of a generic (random) function / to the locations x G X N . Following Equation (jT|), the max- 
stable random vector Z(x) is generated by the Poisson point process = {Qipi(x),i > 1} 
whose intensity measure on the cone E = [0, oo) N is denoted by //,.. Due to its homogeneity, 
the exponent measure /j, x can be factorized into a radial part on (0, oo) and an angular part on the 

unit Li-sphere Sn-i = {z G E : \\z\\ = 1}, where ||z|| = z\ -f-1- zn, for z = (z l7 ..., zn) G 

E. More precisely, a change to polar coordinates under the map U : E —> (0, oo) x SV-i, 
U(z) = (IMU/IMI) yields 


Px(F)= / Vx°U 1 (dr, ds) = N / r 2 drH(ds ), 

JU(F) Ju(F) 


(7) 


for any Borel subset F C E. The probability measure H on Sn -i is called spectral measure of 
Z(x) and it satisfies 


= JV ‘ 1 > J = 1.JV- 

Equation (|7]) shows that we can represent the process d> x as 

= {U^ 1 (R i ,Q i ) :*>!} = {RiQi ■ i > 1}, 
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where {R L : i > 1} is a Poisson point process on (0, oo) with intensity Nr~ 2 dr and Qi , 
i > 1, are independent samples from the spectral measure II on SV_ i. The great advantage of 
this representation is that the components of Qi are bounded by 1. This ensures that Z(x) = 
maxj>! RjQ, can be simulated exactly by generating the largest R, first until no more of the 
remaining points RiQi can contribute to the maximum. 

The only difficulty is thus to generate the random variables Qi from the probability measure H 
on the (N — 1)-dimensional positive sphere Sn-i- The following theorem gives the solution to 
this problem for the max-stable distribution Z(x). 


Theorem 1. Let T r , i > 1, be independent copies of a random variable T with uniform dis¬ 
tribution on the discrete set {1,..., N}. Further, for any k = 1,..., N, let Y^ k \ i > 1, be 
independent random processes with distribution P Xk as in (|4]). Then, the SN-i-valued random 
variables 


Qi 


YpHx) 

ik <t,i mii’ 


i > 1, 


are independent with distribution H. Consequently, with { R,, i > 1} as above, 

Y^ Ti) (x) 

Z(x ) = max R .,—-. 

ei II if" Mil 

Proof. For any k — 1,..., N, Eq. ([4]) implies 

fc f( X k)l{f(x)/\\f(x)\\eA} ^(d/) = f c ^{f(x)/\\f(x)\\£A} Pxki^f). (8) 

We compute the p , x -measure of the set f/ _1 ((u, oo) x A) for u > 0 and a Borel set A C SV-i ■ 

P POO 

dx(U ((«,oo) x A)) ^{CII/f^lllx^l^f/f^l/ll/f^lllsAlC dC u(df) 

Jc Jo 

1 i N 

= ~ / ll/(^)ll 1 {/W/||/( a; )||eA}^(d/) = /(• rfc ) 1 f/( a: )/ll/( a: )lleA} v(df) 

u Jc u , Jc 


k =1 JC 
N 


^' r N 1 ^ r 

/ l {f{AI\\RA\\&A} Px k (df) = ~ ' JfYl / hfW/WfWUcA} Px k (df), 

k =1 JC k =1 JC 


where the second last equation follows from ([8). Let Y {k \ k = 1..... A r , be independent ran¬ 
dom processes with distribution P Xk , respectively, and let T be an independent uniform random 
variable on {1,..., N}, then the above implies that 

p x (U-\(u, oo) x A)) = ^ • pr{ ||y (r) ^|| e a| . 

Comparing this with ([7]) yields the assertion of the theorem. □ 


Theorem [I] shows how to simulate from the spectral measure H. It requires only to be able to 
simulate from the distributions P Xk , k = 1,... ,N. Algorithm [Tj an adaptation of Schlather’s 


[2002] algorithm, provides an exact sample from the max-stable process Z at locations x. 
REMARK 2. The results on the distribution P Xk for stationary Brown-Resnick processes ob¬ 


tained in Subsection 5.2 reveal that Algorithm^ is identical to the algorithm by\Dieker and 
Mikosch [\20I5\I in this case. 
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Algorithm 1. Simulation of a max-stable process Z, exactly at x = (xi,..., x N ) 

1 Simulate <C _1 ~ Exp(lV) and set Z(x) = 0. 

2 While (£ > min(Z(xi),..., Z(xn))) { 

3 Simulate T uniform on {1,..., N} and Y according to the law P XT . 

4 Update Z(x) by the componentwise max(Z(x), CE(x)/||U(x)||). 

5 Simulate E ~ Exp(iV) and update by + E. 

<* } 

7 Return Z. 


3.3 Simulation via extremal functions 

We now introduce the second procedure for exact simulation of the max-stable process Z at 
locations x = (xi,..., xjv) £ X N . For n = 1,..., N we consider the extremal and subextremal 
point processes <f>+ = )Xn} and We have <f>+ = {0+}i <i< n but there 

may be some repetitions in the right-hand side. We define the 72th-step maximum process 

Z n {x) = max 0(x) = max 0+(x), x e X. (9) 

l<i<n 

By the definition of extremal functions we have Z{xf) = 0+(x0 and clearly 

Z(xi) = Z n (xi), 7 = 1,..., n. (10) 

Hence, in order to exactly simulate Z at locations x, it is enough to exactly simulate <3>^. We will 
proceed inductively and simulate the sequence (of )i< n <jv according to the following theorem. 

Theorem 2. The distribution o/’(0+ji< n </y is given as follows: 

• Initial distribution: the extremal function 0+ has the same distribution as F\ V) where F\ 
is a unit Frechet random variable and Y\ an independent random process with distribution 
P Xl given by ©• 

• Conditional distribution: for 1 < n < N — 1, the conditional distribution of(f)f +i with 
respect to (4> x )\<i<n is equal to the distribution of 


^n+1 


argmax^U.^ ( t ) ( x n+ 1) tf®n +i f 0 

argmax 0g$ + 0(x n+ i) if$ n + i = 0 


where $ n+ i is a Poisson point process with intensity 


l<i<n}^-{f(x„+i)>Z n (x n+1 )}Md/) 


( 11 ) 


and Z n is defined by ([9]). 

Proof The distribution of 0+ is given in Proposition |TJ We prove the result for the conditional 
distribution of 4>f n+1 with respect to (0+)i<j< n . Recall that $0 = {0+,..., 0+ } and that 
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according to Proposition [3j the conditional distribution of with respect to <f>+ is equal to the 
distribution of a Poisson point process with intensity 

l<i<n}/^(d/) — 1 {f(xi)<Z n (xi), l<i<n}A i (d/), (12) 

where the equality follows from Eq. (jTO]). In order to determine 4>t n+1 we focus on the functions 
0 G satisfying <j)(x n+ 1 ) > Z n (x n+1 ) and consider the restriction 

^n+1 n {/ e C, f(x n+ 1 ) > z n (x n+1 )j. 

It follows from Eq. that conditionally on (^>+)i<j< n , $ n+ i is a Poisson point process with 
intensity given by Eq. ( [TT| ). We distinguish two cases: 

• if $ n+ i = 0 then there is no function in <3>“ exceeding Z n at point x n+ i, that is, Z(x n +i) = 
Z n (x n + 1 ) and 0+ +1 = argmax^+ </>(x n+1 ). 

• If <!>„+] 0 0 then there is some function in <3>“ exceeding Z„ at point x n+1 , that is, 
Z(x n+ 1 ) > Z n (x n+ 1 ) and 0+ +i = argrnax^^ 0(x n+ i). 

This concludes the proof of Theorem [2} □ 

From the above theorem, one can deduce Algorithm [2] for exact simulation of the max-stable 
process Z at locations x = (aq,..., x N ). According to Proposition [5] and Remark [lj the dis¬ 
tribution P Xn+ 1 can be used to simulate < 3> n+ i with intensity ( [TT| ). Hence, as for the spectral 
method, the second algorithm requires only to be able to simulate from the distributions P Xn , 
n — 1,..., N. Figure [2] illustrates the procedure. 


Algorithm 2. Simulation of a max-stable process Z, exactly at x = (aq,..., x N ) 

1 Simulate C -1 ~ Exp(l) and Y ~ P X1 . 

2 Set Z(x) = C Y(x). 

3 For n = 2,..., N: 

4 Simulate <C _1 ~ Exp(l). 

5 while (C > Z(x n )) { 

6 Simulate Y rsj P Xn . 

7 If (Y(xi) < Z(xi) for alii = 1, • • • , n — 1, 

8 update Z(x) by the componentwise max(Z(a:), (Y(x)) . 

9 Simulate E ~ Exp(l) and update £ _1 by C" 1 + E. 

10 } 

li Return Z. 


4 Complexity of the Algorithms 

In this section, we aim at assessing the complexity of Algorithms |T] and [2] as a function of the 
number N of simulation sites. Both algorithms contain the simulation of exponential random 
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Figure 2: Simulation of Z via Algorithm [ 2 ] at locations (xi, x 2 , x 3 , x 4 ). Initial process 0+ is 
always accepted (top-left). Second process 0+ is accepted as it exceeds 2j = 0+ at x 2 but not 
at ,x'i (top-right). Third process 0+ is equal to 0+ since $3 = 0 (bottom-left). First sample of 
P X4 (grey line) is rejected since it exceeds Z 3 at x 3 ; second sample is valid and thus called 0+ 
(bottom-right). 


variables E and the simulation of A r -dimcnsional random vectors Y (x) according to a mixture 
of the laws P X1 ,..., P XN . The simulation of E causes much less computational effort than the 
simulation of Y and can therefore be neglected in the analysis of the algorithmic complexity. We 
thus consider the number C\(N) and C 2 (N) of random vectors Y(x) that need to be simulated 
by Algorithm [I] and [^respectively to obtain one exact simulation of Z{x). Interestingly, one can 
provide simple expressions for E(Ci(iV)) and E(G' 2 (A r )). 


Proposition 4 . The expected number of random vectors Y(x) that are needed for exact simula¬ 
tion of Z at x = (xi,. .., xn) are: 

Algorithm 1: E(Ci(AT)) = A^E ^nrax Z(xj) -1 j 

Algorithm 2: E(C' 2 (AT)) = N 

Furthermore, E(C'i(A/’)) > E(C , 2 (A/’)) with equality if and only if Z{xf) = . .. = Z(xn) almost 
surely. 


The expectation of CfN) can be calculated similarly to Proposition 4.8 in Oesting et al. [20131. 


9 














The proof for the expectation of C 2 ( N) is more difficult and, for the sake of brevity, it is provided 
as a supplementary material to this paper. 

We conclude this section with some comments on the complexity of our algorithms and a com¬ 
parison with other exact simulation procedures. Proposition [4] shows that, for any max-stable 
process, Algorithm [2] is more efficient than Algorithm [T] in terms of the expected number of 
simulated functions. As the spectral functions follow either one of the laws P Xl ,..., P XN or a 
mixture of these, the simulation of a single spectral function is equally complex in both cases. 
Thus, the new Algorithm [2] is always preferable to the generalized Dieker-Mikosch algorithm, 
Algorithm |Tj 

Next, we compare Algorithm [2] with the exact simulation algorithm via the normalized spectral 
representation proposed by |Oesting et al.| [ ]2013p . By Proposition 4.8 in |Oesting et al.| [ |2013[ , 
the number C-s(N) of simulated normalized spectral functions satisfies 


E(C 3 (N)) — ( [ max'0( a T) z/ (d'0) J E (maxZ^j 


2=1 


\-l 


\ 2=1 


and, thus, depends both on the geometry of the set {xi ,..., x n} and on the law of the max- 
stable process Z. For simulation on a large and dense subset of X, the algorithm via the nor¬ 
malized spectral representation is more efficient than Algorithm [5] as E(C 3 (iV)) is bounded by 
(/ sup xeA . ip(x)iy(d‘ip)) E (sup xgA - Z(x ) _1 ) while EC 2 (iV) = N grows with the size of the sub¬ 
set. If N is small or moderate and the vector Z(x) is weakly dependent, we may also have 
E(C 2 (N)) < E (C 3 (N)). 

Besides the efficiency in terms of the expected number of simulated functions, we also need 
to take into account the complexity of the simulation of a single spectral function. Exact and 
efficient simulation procedures for the normalized spectral function are known for some cases 
only (such as moving maxima processes), while they are not available in other cases like Brown- 
Resnick or extremal-t processes. In contrast, the random functions in Algorithms [T] and [2] with 
distributions P xo in Q, x 0 G X, can be simulated efficiently for the most popular max-stable 
process and extreme value copula models. Indeed, in Section [5] below we provide closed form 
expressions for various important examples. 


5 Examples 

5.1 Moving maximum process 

The parameter space is X = 7L d or W 1 and A denotes the counting measure or the Lebesgue 
measure, respectively. A moving maximum process on A is a max-stable process of the form 


Z(x) = max Qh(x — Xi), x G X, 

i> 1 


(13) 


where {(0, x*), i > 1} is a Poisson point process on (0, oo) x X with intensity measure C -2 dC x 
A(dy) and h : X —> [0, oo) is a continuous function satisfying f x h(x)X(dx) — 1. A famous 
example is the Gaussian extreme value process where h is a multivariate Gaussian density on 
R d flSmithL|1990|. 
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Proposition 5. Consider the moving maximum process ( [T3] ). For all xq e X, the distribution 
P XQ is equal to the distribution of the random function 


hf + x -x 0 ) 

Kx) 


with x h(u)X(du). 


All proofs of this section can be found in the supplementary material to this paper. 


5.2 Brown-Resnick process 

We consider max-stable processes obtained by representation (jT[) where v is a probability mea¬ 
sure on C given by 


u{A) = pr [exp (Wf) - a 2 f)/ 2) E A] , A C C Borel 


(14) 


with W a sample-continuous centered Gaussian process on X with variance a 2 (x) = E [W{x) 2 }. 
In other words, v is the distribution of the log-normal process Y(x) = exp ( W(x) — cr 2 (x)/2), 
x E X. The relation E[exp{lH(x)}] = exp(cr 2 (x)/2) ensures that E[Y(x)] = 1 and Equation 
([2]) is satisfied. 

An interesting phenomenon arises when X = 7L d or and W has stationary increments: 


Kabluchko et al. [2009] show that the associated max-stable process Z is then stationary with 


distribution depending only on the semi-variogram 


~,(h) = -E[{W(h)~W(0)Y}, hex. 

The stationary max-stable process Z is called a Brown-Resnick process. However, our results 
apply both in the stationary and non-stationary case [cf., Kabluchko 2011] and except stated 
otherwise we do not assume that W has stationary increments. 

Proposition 6. Consider the Brown-Resnick type model a For all .i'o E X, the distribution 
P xo is equal to the distribution of the log-normal process 


Y(x) = exp ( W(x) - W{x 0 ) - -Var[lE (x) - W(x 0 )] 


x E X. 


Remark 3. It is easy to deduce from the proposition that in the Brown-Resnick case where W 
has stationary increments, then Y has the same distribution as 


exp I IE (x — xo) — W (0) — ~ x o, x E X. 

Remark 4. The finite dimensional margins of Brown-Resnick processes are Hiisler-Reiss dis¬ 
tributions [cf., Hiisler and Reiss 1989 ] and the above therefore provides a method for their exact 
simulation. 
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5.3 Extremal-t process 


We consider the so called extremal-/ max-stable process [cf., Opitz 20131 defined by represen¬ 
tation ([]} with v the distribution of the random process 


Y(x ) = c a max(0, W(x)) a , x G X, 


(15) 


where a > 0, c a = Ti 1 / 2 2~^ a ~ 2 ^ 2 /T (-Mp), and W a sample-continuous centered Gaussian pro¬ 
cess on X with unit variance and covariance function c. The constant c a is such that E[F(x)] = 1 
so that Equation (|2|) is satisfied. For a = 1, the corresponding max-stable process in (]Tj) coin¬ 
cides with the widely used extremal Gaussian process by |Schlather[ [|20021. 


Proposition 7. Consider the extremal-t model G3- For all Xq G X, the distribution P X() is equal 
to the distribution o/max(T, 0)", where T = (T{x)) x& x ls a Student process with a + 1 degrees 
of freedom, location and scale functions given respectively by 


/i(x) = c(x o,x) and c{x i^xf) 


c(xi,x 2 ) - c(xq,xi)c(xq, xf) 
(cr + 1) 


5.4 Multivariate extreme value distributions 


In this section, we review some popular models for multivariate extreme value distributions, i.e., 
the case when X = {1, ..., N} in ([[]) is a finite set for some fixed A r G N. For these models, we 
explicitly calculate the measure Pj 0 for any j 0 = 1,... ,N. Unless otherwise stated, all random 
vectors are iV-dimensional in this section. Multivariate extreme value distributions differ from 
extreme value copulas only by a change in the marginal distribution, so that our methodology 
applies directly to exact simulation of extreme value copulas. For more details on the models, 
we refer to Gudendorf and Segers p010 |. 


Logistic model 

The symmetric logistic model in dimension N with parameter 9 G (0,1) corresponds to the 
max-stable random vector with cumulative distribution function 


pr[Z < z] = exp (X ] i=1 z j 1/6 ) J ’ 2 = (zi,... ,z N ) G ( 0 , 00 )^. (16) 

Proposition 8. Let /3 = 1/6. In the logistic model ( p~6| ), the probability measure Pj 0 for any 
jo — 1, ■ ■ ■, N is equal to the distribution of the random vector 




where F\, .... Fn are independent, Fj, j f jo, follows a Frechet(/5, cf) distribution with scale 
parameter cp = T(1 — 1 /and (Fj 0 /cf)~ 13 follows a Gamma(l — 1//3,1) distribution. 


REMARK 5. The asymmetric logistic distribution can be represented as the mixture of symmetric 
logistic distributions; see Theorem 1 in Stephenson / 2003]j , for instance. As a consequence, 
Proposit ion [5] also enables exact simulation of asymmetric logistic distributions. 
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Negative logistic model 

The negative logistic model in dimension N with parameter 0 > 0 corresponds to the max-stable 
random vector Z with cumulative distribution function 


'jeJ 


z. 


- 1/8 


, zG(0, oo) 


N 


pr[Z < z] = exp £ (£ 

\0^JC{l,...,iV} / 

Proposition 9. In the negative logistic model 0. the probability measure Pj 0 for any j 0 
1 ,,N is equal to the distribution of the random vector 

Wi W N 


(17) 




\- ■ J U 'JO / 

where W\,..., Wn are independent, Wj, j f jo, follows a Weibull(0, cf) distribution with scale 
parameter cq = T(1 + 1 /6)~ l and (Wj 0 /co) e follows a T(1 + 1/9, 1) distribution. 

Dirichlet mixture model 

The Dirichlet mixture model was introduced by Boldi and Davison [20071. In dimension N, the 
model corresponds to the max-stable random vector given by 

Z = ma xQ(NYi) (18) 

i> 1 

where the Yf s are independent identically distributed random vectors on the simplex 

S N -i = [y e [0, l] n : J2 j=1 yj = 1 }- 


The 

form 


l. -- .7 = 1 J 

distribution of each Yi is a mixture of m Dirichlet models, i.e. its Lebesgue density is of the 

l 

m 

h{y ) = ^ 7r fc diri (y \ a lk ,..., a Nk ), y = (w u ..., w N ) e SV-l, (19) 

k =1 

where 7T k > 0, k — 1,..., m such that n k — 1, &ik > 0, i — 1,..., N, k — 1,..., m, and 

i A „,_i , nf=i r («i) 


d Hv \<*i ,...,a N ) = y* j \ B(a) = —^ 

Here, the parameters n k and a ik , % — 1,..., N, k — 1,..., m, are such that 

= XT’ j = 


N 

7=1 a P 


( 20 ) 


m 

E K1 = £ ^ 


k =1 


a jk 
L^i =1 


Proposition 10. In the Dirichlet model ©. we have for any j 0 = 1 ,,N that Pj 0 = 
J2kLi ^kP-jo* where hk = Ttk<^j 0 k/{Y^iL i anc ( Pjf* equalto the distribution of the random 


vector 


G \ 


<k) 


G 


-(*)' 

N 


dk) ’ • • • ’ r (k) 

X r jo ^jo / 

(A;) (k) 

and G\ ,..., G N are independent random variables with Gamma distribution 


Gj^ ~ Gamma(o!j 0 fc + 1,1) and Gj ~ Gamma(a^, 1), j jo. 
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6 Simulation on dense grids 

In many applications, one is interested in simulating a max-stable process Zona dense grid, 
e.g. x = X fl ( eZ) d . As discussed in Section |4[ on average, this requires the simulation of 
EC^AT) = N random functions in Algorithm [2| that is, the simulation of N random vectors 
of size N. For small e, N will be large and the procedure can become very time consuming. 
Thus, one might be interested in aborting Algorithm [2] after m < N steps, ensuring exactness of 
the simulation only at locations Xi,..., x m . In this case, an alternative design of the algorithm 
which efficiently chooses the subset of m locations might improve the probability of an exact 
sample at all N locations. 

For comparison of two designs, we introduce the random number 

Nq = min{m e {1,... ,N} : Z m (x) = Z N (x)}. 

For n > N 0 , the algorithm does not provide any new extremal functions, but all the simulated 
functions are rejected. Hence, N 0 is the optimal number of iterations before aborting the algo¬ 
rithm. One design is preferable to another if its corresponding random number N 0 tends to be 
smaller. An efficient design should thus simulate the extremal functions at an early stage of the 
algorithm. Based on the intuition that 0+ is likely not to be contained in if Z n {x n +\ ) is 
small, we propose the following adaptive numbering xW,..., x (N> of points in Algorithm]^ 

set x^ 1 ' = X\ and x *- n+1) = argmin { Z n (x ) : x E {xi,..., Xn} \ ..., x^}} , (21) 

for n = 1,..., N — 1. A simulation study indicates that this adaptive version is clearly preferable 
to Algorithm [2] with a deterministic numbering of locations. The advantage is particularly big in 
the case of strong dependence which corresponds to simulation on dense grids. More details on 
the simulation study and its results are provided in the supplementary material to this paper. 
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A Supplementary material 
A.l Proof of Proposition [4] 

Proof of Proposition^ In order to analyze the complexity of Algorithm [2j we consider each 
step of the algorithm separately. In the nth step, i.e. for sampling the process perfectly at site 
x n , we simulate Poisson points £ and stochastic processes Y, until one of the following two 
conditions is satisfied: 

(a) £ < Z n - \ (x n ). This condition is checked directly after the simulation of £ and, in this 
case, no stochastic process Y needs to be simulated. 

(b) £ > Z n _i(x n ) and £ Y{xf) < Z(xi) for all 1 < i < n — 1. In this case, Z is updated and 
£Y is an extremal function as it contributes to Z at site x n (and possibly also at some of 
the sites x n+ i ,..., x at). 


Thus, any stochastic process that is simulated is either rejected, i.e. it is not considered as con¬ 
tribution to Z as it does not respect all the values Z(xf ),..., Z(x n _i), or it leads to an extremal 
function. Denoting by > 1} a Poisson point process on (0, oo ) x C with inten¬ 

sity measure £ _2 d£ P Xn ( df>), the random number C 2 (N) of processes simulated in Algorithm [2] 
satisfies 


N 


£7 2 (JV) = |*+.„,! + £ 


n =2 


*>1: d n) > Z(x n ), £ 


(n) 


n ~ 1 Z(x 
> min , , 1 

i - 1 V>”’( 



( 22 ) 


In this formula, the term 1<3>f , I is the number of extremal functions that need to be sim- 

X'E 1 J • • • j J 

ulated, and the term with index n in the sum is the number of functions that are simulated but 
rejected since > Z(xf) for some j <n — 1. For the computation of the expectation 

of the second term, conditionally on ,, i.e. for fixed Z(x 3 ), 1 < j < n — 1, the two 

sets 


$i n) = {(d n \ i n) ) ■ d n VI n) (^) > Z(xj) for some j = 1, ..., n - 1} 
and ^ : ^ n) ^ n) (xj) < Z( Xj ) for all j = 1,..., n - 1} 

are independent Poisson point processes with intensities £ _2 l| ?>min „-i^ ( . a , ,y^ x .))}d£ P Xn (dip) 

and £ _2 1 j €<min nji (z(a . j)/v , (3;i)) t.d£ P Xn (dip), respectively. Conditioning further on d4 n) , Z(x n ) is 
also fixed and we obtain 


E 

= E (E 
= E 


(d” , A" ) ) : d”’ > %), d”' > min 


(n) ^ ™ •£ Z(Xj) 


j ~ 1 4 n) ( x j) 


(£,^)G$! n) : £>Z(x n ) 


$ 




*?’)) 





£>min 


1 P Xn (df) 

i>iff) j 
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= E 


^min 


1 n- 1 Y n (Xj)\\ 

ZiXn)^ Z( Xj ) j) 


where Y n ~ P Xn and Z are independent. The relation min {a. b} = a + b — max{a, b}, a, b E M, 
and the fact that Y n {x n ) = 1 almost surely yield 


E 

= E 


/ (-(n) / (n)\ n ) ry t \ c (n) n Z Z^Xj) 

{Q 0 : Q > Z{x n ), Q > mm / 

J 1 W A x i) 

i n— i Y n {xj)\ ( n 

+E| % a . x zfeyJ- E l, i ;s £ zfe)J 

“ E l*ta ..,)!■ 


Z(x r 

= l + E|$+ iv ^_ i} 


as E|<h|' a , i t , | = E (max" =1 Y n {xf)/Z{xj)) by Lemma 4.7 in Oesting et al. 2013]. Thus, by 


fl22|), we obtain 


E C 2 {N) = E|$+ 






= JV-l + E|*+ l) | = JV. 


Moreover, by ([2]), we have that E Z[x,i) 1 = 1 for i — 1,..., N, and, thus, 

E ^niax Z(xi)^ 1 j > 1, 


with equality if only if Z(x i) = ... = Z(x^) holds almost surely. 


□ 


A.2 Proofs for Section [5] 

A.2.1 Moving maximum process 

Proof of Proposition^ In the case of the moving maximum process ( fT3j ), the measure v asso¬ 
ciated with the representation (|T|) is 


z'(A) = / l {M ._ x)eA }A(dy), A c C Borel. 

J x 

We deduce from Proposition^ 

p x 0 (A) = l {f/f{xo)eA} f{x 0 )u(df)= / l {H ._ x)/hixo - x )eA}h(xo - x)A(dy) 


' x 


lx 


Mt— xo)/h(u)£A} h[u) A(dtt) 


where the last line follows from the simple change of variable xq~x = u - This proves the result 
since h{u)\(du) is a density function on X. □ 


18 














A.2.2 Brown-Resnick process 

Our proof of Proposition [6] relies on the following lemma on exponential changes of measures 
for Gaussian processes. 

Lemma 1. The distribution of the random process {W (x)) xe x under the transformed probabil¬ 
ity measure pr = e w(x T- rj (U)/^ f jp r Is equal to the distribution of the Gaussian random process 


W(x) + c(xo,x), x G X, 


where c(x, y ) denotes the covariance between W(x) and W(y). 

Proof of Lemma\Jj We need to consider finite dimensional distributions only and we compute 
for some x ±,..., x k G X the Laplace transform of (W (x;))i<i<fc under the transformed proba¬ 
bility measure pr. For all 6 = {9 1 ,..., Of) G we have 


c(e u ...,o k ) = E 




= E 


,ir(io)-o' 2 (io)/2„Et 1 OiW(xi) 


= exp ( ^9 T J39 - ^a 2 (x 0 ) ) , 


(23) 


with 9 = (1, 9) G R k+1 and E 
block decomposition 


(c(xi, Xj)) 0<i j <k the covariance matrix. We introduce the 

g = ( v 2 (x 0 ) E o, fc \ 

V S / 


with E = 


c[x. 






and E 


k, 0 


_ V T 

— ^0 .k 


(23 ) can be rewritten as 


(c(x 0 , Xi))i<i<k- The quadratic form in Equation 


\~9 T t~9 - lo 2 (x 0 ) = \ (o 2 (x o) + 9 T Z9 + 26 ,T E/ Cj0 ) - la 2 (x 0 ) = 9 T E*, 0 + Vs0. 


We recognize the Laplace transform of a Gaussian random vector with mean E/,. :0 and covariance 
matrix E whence the Lemma follows. □ 


Proof of Proposition^ Equations 0 and ( |T4j ) together with Lemma [T] yield, for all Borel set 

AcC, 


P X0 (A) = / l {f/f ( x ) eA} f(x)v(df) = E 


e W(x 0 )-^cr 2 (xo)i 


| e W(-)-^cr 2 (-) y e ^( a: o)-2' :r2 ( x o)g J 4| 


= pr 

= pr 

= pr 


exp ( W(-) - W(x 0 ) - -(a 2 (-) - a 2 (x 0 )) ) G A 


exp ( Wf) + c(x 0 , •) - W{x q) - c(x 0 ,a:o) - -(<r 2 (-) - a 2 (a:o)) ) G A 


exp ( W(-) - W{xf) - -(cx 2 (-) + o- 2 (x 0 ) - 2c(x 0 , ■)) ) e 3 


Using the fact that for all x G X 

a 2 [x ) + a 2 (x 0 ) — 2c(xo,a;) = Var[VF(a;) — W(x 0 )] 
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we deduce that P Xo is equal to the distribution of the log-normal process 


Y(x) = exp ( W(x) - W(x o) - -Var [W(x) - W(x 0 )] ) , x G X. 


This proves Proposition |6j 


□ 


A.2.3 Extremal-/ process 

In the sequel, we write shortly z a = max(0, z) a for all real numbers z. 


Lemma 2. The distribution of the random process {W(x) /W{xf) :r(iX under the transformed 
probability measure pr = c a W(x 0 )“dpr is equal to the distribution of a Student process with 
a + 1 degrees of freedom, location parameter p k and scale matrix S^ given by 


Pk ^/s,o and Sk k 


Si. — Sfc.nS 


k,0^0,k 


Oi ~|~ 1 


where S fc = (c(x i ,a; i )) 1 < i>i < fc and S fc>0 = ££* = (c(x 0 , 

Proof of Lemma [2] We consider finite dimensional distributions only. Let k > 1 and xi,..., x k G 
X. We first assume that the covariance matrix S = (c(xi,Xj)) 0<ij<k is non singular so that 
(W(xi)) 0 <i< k has density 


9(y) = { 2vr ) (fc+1)/2 det(S) 1/2 exp ^?/ T S 1 y^j with y = 


Setting z = (yi/yo)i<i<k, we have for all Borel sets A 1} ..., A k c 

'W{xi 


pr 


yw{x 0 ) 


e Ai i = l,... ,k 


lRk+1 


^-{yi/yo&Ai, i=i,...,fc} c cd/o d{y ) d 2/ 


MzitAi, <=i.*} / c Q ?/o^(|/o,2/o-)2/od2/o d - 


We deduce that under pr, the random vector (W(xj)/W(^o))i<t<fc has density 

POO 

9 {z) = / c a yQ +a g(y 0 ,y 0 z) dy 0 

Jo 

_ /*°° 

= c Q (27r)-( fc+ d/ 2 det(S)- 1 / 2 / ^ +Q exp 

Jo 

with 5 = (1, z). Using the change of variable u = zTj i 1 z y% , we get 


zT- 1 ^ 


2/o d ?/o 


?/o +Q exp 


^s- 1 ^ 


2/o d 2/o = x 


1 / fs- 1 ^ 


a+k+1 

2 /»oo 




(k+a- l)/2 exp (_ M ) dM 


1 / FE^z 

2 


a+fc + l 
2 


r 


k CM, H - 1 
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and we obtain after simplification 


r> ( fc+g+1 \ o+fc+i 

g( Z ) = vr - fc /2_i_|_i. det ( S )-l/2 

Introducing the block decomposition E = ( ^ ~f u ' ), the inverse matrix is 

\ ^k,0 Efc 

y~i _ ( 1 + £o,fc(£fc — S fci0 So,fc) _1 ^fc,o — S 0 ,fc(£fc — £fc,oSo,fe) _1 

— (£fc — £fc, 0 £o,fc) - 1 £fc,o (^fe — S fc> o£o,fc ) _1 

By the definition of fi k and E fc , we have 


E _1 = 


1 + a + fj > lH k 1 tx k 
1 + OL \ — E fc 1 /ifc 


.T V — 1 


E* 1 


and 


z T E _1 z = (1, ^) T £ _1 (1, z) = ( 1 + 

Finally, we obtain after simplification 


(z-p. fc ) T E fc jz-fik ) 

OL - 1 - 1 


-t/2 r )—1/2 C J + (* - - H) 


q+fc+1 

2 


g(z) = n t/2 (a + l) k/2 - ' Jf ■' det(S t 


T(2±i) 

We recognize the variate Student density with a + 1 degrees of freedom, location parameter 
/i and scale matrix E/,.. □ 

Proof of Proposition^ Consider the set 

-A = {/ G Co; /(^i) e .Ai, ■ ■ ■ , /(xfc) G Ak}. 

Equations Q and ( [15] ) together with Lemma [2] yield, 

PxM) J ^~{f/f{x)&A} f E \c a W (xo) ^-{W{xi) a /W{xo) a &Ai, i=l....,k}\ 

= pr \W(xi) a lW(xo) a e Ai, i = 1, ... ,k\ 

= pr [T“ G A, i = A:] 


where T = (Ti,..., Tf has a multivariate Student distribution with a + 1 degrees of freedom, 
location parameter p k and dispersion matrix E^. This proves the result. □ 
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A.2.4 Multivariate extreme value distributions 


Logistic model 

Proof of Proposition^ It is easily shown that the logistic model admits the representation 

Z = max QF,j 

i> 1 


where the Ff s are independent random vectors with independent Frechet(/3, cf )-distributed 
components. To check this, we compute 


E 


N Fj 
max — 
3 =i 


i o 


IF, 


pr 

-( ZjU/ci 


N Fj 

max — > u 


3 =i 2 ; 


3 


d u = 


r N 


dw = 


3 -^£f=l {Zi/ep) 


n =1 < z 3 u \) du 

d u = 




Z-^j=l 3 ) 


For the computation of the last integral, we recognize the expectation of a Frechet distribution. 
Next we use the fact that P ]0 is the distribution of F/Fj 0 under the transformed density 



-1-/3 

e -(yk/c{)) p ' 


We recognize a product measure where the jth margin, j f j 0 , has a Frechet (/T cf) distribution. 
The joth marginal has density 

C P\ C PJ 

and a simple change of variable reveals that this is the density of with Z ~ Gamma(l — 

1//5,1)- □ 

Negative logistic model 

Proof of Proposition^ Similarly to the logistic model, we have the spectral representation 

Z = max(j Wj 

i> 1 


where the IF,;’s are independent random vectors with independent Weibull^, c^)-distributed 
components with scale parameter cq = | M 1 | /0j . To check this, we compute 


r n Wj 1 

POO 

\ N Wj ) 

max —- > u 

max—- 

= / P r 

. J =1 Z 3 . 

Jo 

. i =1 Z J 


N 

1 ~ n ( X _ e- ( ^“ /cs)e ) 

3=1 



poo / N 

du = / II — TTpr[fFj 
J o V 3=1 

)*.= 


< Zju\ ) du 




- I/ V(l + l/9) = - ^(-1)W (£, ' A ~ I/ " 

J 
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For the computation of the last integral, we recognize the expectation of a Weibull distribution. 
As for the logistic model, Pj 0 is the distribution of W/Wj 0 under the transformed density 



g-(yfc/ce)®_ 


We recognize a product measure where the jth margin, j f j 0 has a Weibull(0, cq) distribution. 
The joth marginal has density 


Vjo 


ce 


Vjo 

ce 


e-i 


~(yj Q / c e) e 


and a simple change of variable reveals that this is the density of cgZ l ^ e with Z ~ Gamma(l + 

1 / 0 , 1 ). □ 


Dirichlet mixture model 


Proof of Propositioning By definition, //„ has the form 


m n. 

Pjo ( A ) = N Y] / yj 0 1 { y /y jo e a } diri (y \ a lk , ..., a Nk ) d y 

i i 0 Sn-i 

■fsjv-i yjo^-{y/yj 0 &A}difi(y \ ctjfc,..., Qfjvfc) dr/ 


k =i 
m 


n Y. 


, A c (0, oo)^. 


k =1 


fs N _, % 0 diri(r/ | ai fc ,..., a Nk ) dr/ 


Thus, //„ is given as the mixture P J0 = J///, n k P^\ where for each k = 1,..., m, the prob¬ 
ability measure P A> is equal to the distribution of the random vector /Y^\ and Y ik) has 

a transformed density proportional to r/j 0 JX/li vf' ~ ' • We recognize the Dirichlet distribution 
with parameters ai k ,..., aNk given by 


— ^jok d - 1 and oij k — oij k j jo¬ 


lt is well known that Dirichlet distributions can be expressed in terms of Gamma distributions. 
More precisely, we have the stochastic representation 


N 


Y ik) = G[ k) 



V i=i j =i / 

where G lk) are independent Gamma(ay fc , 1) random variables. The result follows since P^ is 
the distribution of Y {k '>/Y^ . □ 
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Figure 3: Histogram of N 0 based on 5000 realizations of a Brown-Resnick process associated 
to the semi-variogram 7 (h) = c||/i|| a with c = 1 and a = 1.5 (left), c = 2.5 and a = 1 (middle) 
and c = 5 and a = 0.5 simulated via Algorithm [2] with the deterministic design (grey) and the 
adaptive design ( [ 21 ] ) (white), respectively. 


A.3 Simulation Study 

We perform a simulation study to compare the adaptive version of Algorithm[ 2 ]introduced in ( |2Tj ) 
to a version, where the numbering of locations is deterministic. The simulation study is based 
on 5000 simulations of a Brown-Resnick process associated to a semi-variogram of the type 
7 (h) = c\\h\\ a on the two-dimensional grid (0.05, 0.15,..., 0.95} x {0.05, 0.15,..., 0.95}. We 
run Algorithm [2] with the deterministic design (the grid points are ordered by their coordinates 
in the lexicographical sense) and with the adaptive design ( [2~i~[ ). The simulation is repeated 
for different values of the parameter vector (c, a) representing strong dependence ((c, a) = 
(1,1.5)), moderate dependence ((c, a) = (2.5,1)) and weak dependence ((c, a) = (5,0.5)). 
The histograms of N 0 are shown in Figure A.3[ For each of the three parameter vectors, the 
number N 0 for the adaptive design is stochastically smaller than the corresponding number for 
the deterministic design. 
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